Rt per provincia¶

Rt di COVID-19 stimato per le province italiane.

Max Pierini, Alessio Pamovio & NOTIZIÆ Telegram Channel

NB: gli $R_t$ stimati da EpiDataItalia sono parametrizzati attualmente sui soli nuovi casi riportati al Dipartimento di Protezione Civile e non sono sovrapponibili al numero di riproduzione stimato da ISS-EpiCentro che utilizza un modello EpiEstim (A new framework and software to estimate time-varying reproduction numbers during epidemics, Cori-Ferguson-Fraser 2013) per il quale sono necessari dati attualmente non disponibili al pubblico (il numero di casi sintomatici con data di inizio sintomi, il numero di casi sintomatici importati da un’altra regione o dall’estero). Per info sul modello in uso da ISS-EpiCentro consultare la pagina FAQ sul calcolo del Rt di ISS. Per chiedere al Governo l'accesso pubblico a tutti i dati grezzi disaggregati di COVID-19 in Italia consigliamo l'adesione alla Petizione DatiBeneComune.

Method¶

A simple method is presented to estimate effective reproduction number $R_t$ of COVID-19 in italian provinces with a Markov chain Monte Carlo and Poisson likelihood parametrized on daily new cases.

Details¶

New cases $y_r$ for each $r$ italian province (source Dipartimento Protezione Civile), will be smoothed with rolling mean (gaussian, window 7, std 2). Smoothed new cases will be adjusted to be $>0$ to avoid negative values (due to data error corrections).

For each day $t$, smoothed new cases $y_{r,t}$ will be supposed distributed as Poisson with $\lambda_{r,t}$ parameter

$$ y_{r,t} \sim \mathcal{P}(\lambda_{r,t}) $$

where $\lambda_{r,t}$ is defined by the serial interval inverse $\gamma$, previous day smoothed new cases $k_{r,t-1}$ and effective reproduction number in time $R_{r,t}$ (ref: Bettencourt & Ribeiro 2008)

$$ \lambda_{r,t} = k_{r,t-1} e^{\gamma (R_{r,t} - 1)} $$

Parameters $R_{{p,t}}$, for each province $p$ and time $t$ (day) will be distributed as half normal with mean equal to previous day posteriors $R_{{p,t-1}}$ and precision $\tau$

$$ R_{{p,t}} \sim \mathcal{{N}}^+(R_{{p,t-1}}, \tau) $$

where, first day $R_{{p,0}}$ (outcome) is set to zero

$$ R_{{p,0}} = 0 $$

and $\tau$ is defined as $1 / \sigma^2$, where $\sigma$ is the standard deviation of $R_t$ priors, previuosly estimated for italian provinces (using the same method but with overarching uninformative prior distribution of $\tau$) as

$$ \sigma \sim \mathcal{N}(\mu_{=0.299}, \sigma_{=0.004}) $$

If previous new cases are zero $k_{p,t-1}=0$, parameter $R_{p,t}$ is undefined, given the chosen function for $\lambda_{p,t}$ parameter of Poisson likelihood, even if it should be $R_{p,t}=0$ (no new cases means null effective reproduction number). Thus, in these cases, priors of $R_{p,t}$ will be forced to

$$ R_{p,t} \sim \mathcal{N}^+(0, \tau) \;,\; k_{p,t-1}=0 $$

and previous new cases will be forced to $k_{p,t-1}=1$, so that $\lambda_{p,t}$ will be

$$ \lambda_{p,t} = e^{\gamma( \mathcal{N}^+(0, \tau) - 1 )} \;,\; k_{p,t-1}=0 $$

model {
    # Overarching Rt standard deviation
    sigma_R ~ dnorm( 0.29929894605552754 , 1 / 0.004218809765237365 )
    tau_R <- 1 / sigma_R^2


    # Serial interval distribution
    SI_mu <- 7.5
    SI_sd <- 3.4
    SI_sh <- SI_mu^2 / SI_sd^2
    SI_ra <- SI_mu / SI_sd^2
    SI ~ dgamma( SI_sh , SI_ra )
    gamma <- 1 / SI

    # First Rt prior
    R[1] <- 0
    for ( t in 2:T ) {
        # Rt prior for k>0
        Rpr[t] ~ dnorm( R[t-1] , tau_R )  T(0,)
        # Rt prior for k=0
        Rnu[t] ~ dnorm( 0 , tau_R )  T(0,)

        # Define Rt prior
        R[t] <- ifelse( k[t-1]==0 , Rnu[t] , Rpr[t] )
        # Avoid k=0 (undefined Rt)
        K[t] <- ifelse( k[t-1]==0, 1 , k[t-1] )

        # Poisson likelihood
        lambda[t] <- K[t] * exp( gamma * ( R[t] - 1 ) )
        y[t] ~ dpois( lambda[t] )
    }
}

model {
    # Overarching Rt standard deviation
    sigma_R ~ dnorm( 0.29929894605552754 , 1 / 0.004218809765237365 )
    tau_R <- 1 / sigma_R^2


    # Serial interval distribution
    SI_mu <- 7.5
    SI_sd <- 3.4
    SI_sh <- SI_mu^2 / SI_sd^2
    SI_ra <- SI_mu / SI_sd^2
    SI ~ dgamma( SI_sh , SI_ra )
    gamma <- 1 / SI

    # First Rt prior
    R[1] <- 0
    for ( t in 2:T ) {
        # Rt prior for k>0
        Rpr[t] ~ dnorm( R[t-1] , tau_R )  T(0,)
        # Rt prior for k=0
        Rnu[t] ~ dnorm( 0 , tau_R )  T(0,)

        # Define Rt prior
        R[t] <- ifelse( k[t-1]==0 , Rnu[t] , Rpr[t] )
        # Avoid k=0 (undefined Rt)
        K[t] <- ifelse( k[t-1]==0, 1 , k[t-1] )

        # Poisson likelihood
        lambda[t] <- K[t] * exp( gamma * ( R[t] - 1 ) )
        y[t] ~ dpois( lambda[t] )
    }
}

DONE!

Time series¶

Latest Rt (IQR)¶

Latest Rt (HPDI)¶

Latest Rt (MAP)¶

-> 2020-11-10 18:34:32.988763
<- 2020-11-10 19:53:20.031080
Completed in 1:18:47.042317